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Abstract. Robust estimation method is proposed to combine multiple observations 
and create consistent, accurate, dense Digital Elevation Models (DEMs) from lunar 
orbital imagery. The NASA Ames Intelligent Robotics Group (IRG) aims to produce 
higher-quality terrain reconstructions of the Moon from Apollo Metric Camera 
(AMC) data than is currently possible. In particular, IRG makes use of a stereo vision 
process, the Ames Stereo Pipeline (ASP), to automatically generate DEMs from 
consecutive AMC image pairs. However, the DEMs currently produced by the ASP 
often contain errors and inconsistencies due to image noise, shadows, etc. The 
proposed method addresses this problem by making use of multiple observations and 
by considering their goodness of fit to improve both the accuracy and robustness of 
the estimate. The stepwise regression method is applied to estimate the relaxed weight 
of each observation. 


1. Introduction 

Since 2007, the NASA Lunar Mapping and Modeling Project (LMMP) has been actively 
developing maps and tools to improve lunar exploration and mission planning [1]. One of 
the requirements for LMMP is to construct geo-registered DEMs from historic imagery. To 
meet this need, IRG has developed the Ames Stereo Pipeline (ASP), a collection of 
cartographic and stereogrammetric tools for automatically producing DEMs from images 
acquired with the Apollo Metric Camera (AMC) during Apollo 15-17 (Figure 1). 

There are many applications of lunar maps, including outreach and education, mission 
planning, and lunar science. The maps and imagery that IRG released in 2009 for “Moon in 
Google Earth” have engendered great public interest in lunar exploration. The cartographic 
products that IRG is currently producing (DEMs, ortho-projected imagery, etc.) will be 
used to plan future missions, to assess landing sites, and to model geophysical processes. 
Consequently, by improving these maps, we directly benefit a wide community. 

The NASA Ames Intelligent Robotics Group (IRG) currently uses a stereo vision 
process to automatically generate DEMs from image pairs [2]. However, two DEMs 
generated from different image pairs have different values for the same point due to noise, 
shadows, etc. in the images (Figure 2a). Consequently, IRG’s current topographical 




(a) The Mapping Cameras System (b) Apollo 15 Orbit 33 

Figure 1: AMC Data System, (a) The AMC captures a series of pictures of the Moon's surface, (b) 
Satellite station positions for Apollo Orbit 33 visualized in Google Moon. 


reconstruction of the Moon contains fairly substantial random errors (Figure 2b). It is 
important to construct consistent DEMs to estimate the photometric properties [3-5]. 

This paper will address this problem by finding robust elevation from multiple DEMs 
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Figure 2: Stereo and Multiple View Correspondences, (a) Two stereo correspondences have different 
elevation values (two crosses), (b) DEMs created by stereo can have substantially large errors. 



that minimize the weighted squared error of all associated DEM patches. The proposed 
method determines the unique 3D position by weighted averaging of all elevation values 
from stereo DEMs. The accuracy and robustness of DEMs produced by IRG will thus be 
improved by making use of multiple observations and by considering their goodness of fit. 
The reconstructed DEMs from lunar orbital imagery are presented. This paper tackles one 
of the challenging problems of planetary mapping from orbital imagery. 

2. Ames Stereo Pipeline 

The Ames Stereo Pipeline (ASP) is the stereogrammetric platform that was designed to 
process stereo imagery captured by NASA spacecraft and produce cartographic products 
since the majority of the AMC images have stereo companions [6]. The entire stereo 
correlation process, from an image pair to DEM, can be viewed as a multistage pipeline 
(Figure 3). At the first step, preprocessing includes the registration to align image pairs and 
filtering to enhance the images for better matching. Triangulation is used at the last step to 
generate a DEM from the correspondences. 

2.1. Disparity Initialization 

Stereo correlation, which is the process at the heart of ASP, computes pixel 
correspondences of the image pair (Figure 4). The map of these correspondences is called a 
disparity map. The best match is determined by applying a cost function that compares the 
two windows in the image pair. The normalized cross correlation is robust to slight lighting 
and contrast variation in between a pair of images [7]. For large images, this is 
computationally very expensive, so the correlation process is split into two stages. (1) The 
disparity map initialization step computes coarse correspondences using a multi- scale 
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Figure 3: Dataflow of the Ames Stereo 
Pipeline. Preprocessing includes the 
registration and filtering of the image pair. 
A stereo correlator (disparity initialization 
and sub-pixel refinement) constructs the 
disparity map based on normalized cross 
correlation. DEMs are generated by a 
triangulation method in which corrected 
camera poses are used by bundle 
adjustment. 


search that is highly optimized for speed (Figure 4c). (2) Correlation itself is carried out by 
sliding a small, square template window from the left image over the specified search 
region of the right image (Figure 4d-f). 

Several optimizations are employed to accelerate disparity map initialization [8]: (1) a 
box-filter-like accumulator that reduces duplicate operations during correlation; (2) a 
coarse-to-fine multi- scale approach where disparities are estimated using low resolution 
images, and then successively refined at higher resolutions; and (3) partitioning of the 
disparity search space into rectangular sub-regions with similar values of disparity 
determined in the previous lower resolution level of the pyramid. 



(a) Left Image (b) Right Image (c) Initialized Disparity Map 



(d) Parabola Refined Map (e) Cauchy Refined Map (f) Bayesian Refined Map 


Figure 4: Stereo Correlation, (a-b) An image pair, (c-f) Horizontal disparity maps, (c) The fast 
discrete correlator constructs the coarse disparity map. (d-f) Refined disparity maps from the 
initialized map. (f) Bayesian sub-pixel correlator generates a smoother map than the others. 



2.2. Sub-pixel Refinement 

Refining the initialized disparity map to sub-pixel accuracy is crucial and necessary for 
processing real-world data sets [9]. The Bayesian expectation maximization (EM) weighted 
affine adaptive window correlator was developed to produce high quality stereo matches 
that exhibit a high degree of immunity to image noise (Figure 4f). The Bayesian EM 
sub-pixel correlator also features a deformable template window that can be rotated, scaled, 
and translated as it zeros in on the correct match. This affine-adaptive window is essential 
for computing accurate matches on crater or canyon walls, and on other areas with 
significant perspective distortion due to foreshortening. A Bayesian model that treats the 
parameters as random variables was developed in an EM framework. This statistical model 
includes a Gaussian mixture component to model image noise that is the basis for the 
robustness of the algorithm. The resulting DEM is obtained by the triangulation method 
(Figure 5). 

2.3. Bundle Adjustment 

After stereo correlation is performed, Bundle Adjustment (BA) corrects the 
three-dimensional postures of cameras and the locations of the objects simultaneously to 
minimize the error between the estimated location of the objects and their actual location in 
the images. Camera position and orientation errors have a direct effect on the accuracy of 
DEMs produced by the ASP. If they are not corrected, these uncertainties will result in 
systematic errors in the overall position and slope of the DEMs (Figure 6a). BA ensures 



(a) Left Image 


(b) Right Image (c) Rendered DEM 


Figure 5: Generation of DEM. (a) and (b) Apollo Metric Camera image pair of Apollo 15 site, (c) A 
DEM of Hadley Rille is rendered from an image pair. 




(a) Initial 3D Mesh 


(b) After Bundle Adjustment 


Figure 6: Bundle Adjustment. Color-mapped, hill-shaded DEM mosaics from Apollo 15 Orbit 33 
imagery illustrate the power of BA. (a) Prior to BA, large discontinuities exist between overlapping 
DEMs. (b) After BA, DEM alignment errors are minimized, and no longer visible. 


that observations in multiple different images of a single ground feature are self-consistent 
(Figure 6b). 

In BA the position and orientation of each camera station are determined jointly with the 
3D position of a set of image tie-points points chosen in the overlapping regions between 
images. This optimization is carried out along with thousands (or more) of similar 
constraints involving many different features observed in other images. Tie-points are 
automatically extracted using the SURF robust feature extraction algorithm [10]. Outliers 
are rejected using the random sample consensus (RANSAC) method [11]. The BA in ASP 
determines the best camera parameters that minimize the reprojection error [12]. The 
optimization of the cost function uses the Levenberg-Marquardt algorithm [13]. 

3. Robust Estimation 


Suppose a set of observations P - {x k } n k=l from a normal distribution with some of them are 
contaminated by outliers. Their Boolean membership to inliers is represented by an 
indicator vector w = [ w fc] e 

0 if x h is outlier 

= h * . (i) 

1 otherwise. 


Let w be the total number of inliers: 


w = £>r ( 2 ) 

1=1 

The membership measure supporting the normality is defined to minimize its square error 
and to maximize the number of inliers. The mean of the inliers is written by 
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The squared error is written by 
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An objective function which minimizes s 2 and maximizes w simultaneously is hard to 
define. Even if it is defined, computing the optimum w is NP-hard. 

The Boolean membership vector is relaxed to be a real value in [0,1] to convert the 
combinatorial optimization problem into a tractable nonlinear optimization problem [14]. 
The membership matrix reflects the likelihood of the inliers. Let us redefine 


w = [ w i\ £ [0,l] n , where w i is the membership of point x i to inliers. Then the equations 
from (3) to (7) are still valid. From the expected value of SSE , 
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The squared error of each point is defined by 
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The statistic to test statistical significance for x t follows the F distribution: 
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The updating principle is to increase the weight of a observation if its p-v alue of the statistic 
is smaller than the significance level or to decrease the weight. Similar with the gradient 
descent method, the update rule is implemented by 


where a and 77 are significance level and learning rate. 

4. Experimental Results 

The NASA Exploration Systems Mission Directorate (ESMD) has been charged with 
producing cartographic products via LMMP for use by mission planners and scientists in 
NASA’s Constellation program. As part of the LMMP, we have produced 70 preliminary 
DEMs and ortho-images derived from Apollo 15 Metric Camera (AMC) orbit 33 imagery 
using the ASP. The DEM Mosaic program is implemented based on the NASA Vision 
Workbench (VW). The NASA VW is a general purpose image processing and computer 
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(a) Robust DEM Mosaic 


(b) Weight Image 


Figure 7: DEM Mosaic and its weight. 


vision library developed by the IRG at the NASA Ames Research Center. 

The unified DEM mosaics are constructed with the stereo DEMs. The significance level 
a was set to be 0.3 and the typical learning rate 77 is 0.01 . Figure 7a shows a DEM mosaic 
image constructed from its adjacent stereo DEMs which has large variances between DEMs 
as shown in Figure 2b. As you can see in the figure, the unified DEM mosaic provides the 
seamless elevation map. Figure 7b shows the estimated weight for the current DEM as we 
can expected that the central part of the DEM has large confidence. 

5. Conclusion 

The robust estimation method was proposed to determine the unique elevation from 
multiple digital elevation models generated by the Ames Stereo Pipeline. Introducing the 
relaxed membership of each observation to be an inlier, the stepwise regression method is 
applied to estimate the relaxed weight of the observation. The proposed method addresses 
this problem by making use of multiple observations and by considering their goodness of 
fit to improve both the accuracy and robustness of the estimate. A gradient descent method 
was used to optimized the weight because it is hard to define the objective function in a 
closed form. The residual analysis will be desirable to provide a quantitative measure of the 
proposed method. The parametric representation of surface model will be valuable to 
enhance the recovery resolution and robustness of the algorithm. 
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